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We introduce a model for the real-time evolution of a rela- 
tivistic fluid of quarks coupled to non-equilibrium dynamics of 
the long wavelength (classical) modes of the chiral condensate. 
We solve the equations of motion numerically in 3+1 space- 
time dimensions. Starting the evolution at high temperature 
in the symmetric phase, we study dynamical trajectories that 
either cross the line of first-order phase transitions or evolve 
through its critical endpoint. For those cases, we predict the 
behavior of the azimuthal momentum asymmetry for high- 
energy heavy-ion collisions at nonzero impact parameter. 



I. INTRODUCTION 

The hydrodynamical model is frequently employed to 
describe multiparticle production processes in hadronic 
collisions [1]. In particular, it predicts characteristic flow- 
signatures as "fingerprints" for non-trivial equations of 
state of hot and dense matter [2]. Such equations of 
state can occur when the effective potential, obtained 
by integrating out some degrees of freedom, exhibits fea- 
tures characteristic of a phase transition in thermody- 
namics [3]. 

More specifically, if there exist two (or more) collective 
states with the same free energy but separated by a bar- 
rier, then behavior characteristic of a first-order phase 
transition may emerge. On the other hand, if no free- 
energy barrier exists, one might expect resemblence to a 
second-order phase transition. This analogy of interact- 
ing quantum field theories with thermodynamics is be- 
lieved to have played an important role in the evolution 
of the early universe [4] and is currently being investi- 
gated in accelerator experiments by colliding beams of 
protons and heavy ions [5]. Classical energy flow and 
hydrodynamic scaling behavior emerges in high-energy 
inclusive processes [1,6] and from the real-time evolution 
of some quantum field theories [7] . 

In this paper, we extend the hydrodynamical transport 
model such that phase transitions related to the restora- 
tion (or breaking) of some global symmetry can be stud- 
ied dynamically. In particular, we focus on chiral symme- 
try breaking at finite temperature [8] , for which we shall 
adopt a relatively simple and tractable phenomenological 
model, i.e. the Gell-Mann and Levy model [9]. 

It has been argued [10-13] that the chiral phase tran- 
sition for two massless quark flavors is second-order at 
baryon- chemical potential its = 0, which then becomes 



a smooth crossover for small quark masses. On the other 
hand, a first-order phase transition is predicted for small 
temperature T and large \ib- If, indeed, there is a smooth 
crossover for fig = and high T, and a first-order transi- 
tion for small T and high /xg, then the first-order phase 
transition line in the (/x^ , T) plane must end in a second- 
order critical point. This point was estimated [10] to be 
at T ~ 100 MeV and fi B ~ 600 MeV (see also [12,13]). 
More recently, it has been attempted to determine the 
endpoint of the line of first-order phase transitions from 
the lattice, using 2+1 quark flavors on N t — 4 lattices [14] 
(see also [15]). Those authors locate the critical point at 
T = 160+3.5 MeV and fi B = 725 ± 35 MeV. Note, how- 
ever, that a reliable extrapolation to the continuum limit 
and to physical pion mass has not been attempted so far. 

There is an ongoing experimental effort to detect that 
chiral critical point in heavy-ion collisions at high ener- 
gies. Note that both, high temperature and high baryon 
density are required to have dynamical trajectories in 
heavy-ion collisions pass reasonably close by the critical 
point. Some dynamical computations [16] of the energy 
deposition and baryon stopping process during the initial 
stage of head-on collisions of large nuclei within semi- 
realistic multi-fluid dynamical models suggest that the 
required conditions may be reached in the central region 
of collisions at Ei a \, ~ 20 — 80 AGeV on a fixed target, or 
in the fragmentation regions of collisions at higher ener- 
gies. However, to our knowledge there has been so far no 
attempt to describe hydrodynamical expansion of the hot 
and dense droplet produced initially for dynamical tra- 
jectories close to the critical point. This paper represents 
an attempt in that direction. 



II. THE MODEL 

In this section we shall present our model for the dy- 
namics of a droplet of quarks and antiquarks, starting 
at high temperature in a state with (approximately) re- 
stored chiral symmetry, and evolving towards a state 
where the symmetry is spontaneously broken. The 
quarks will be described as a heat bath in local ther- 
mal equilibrium that evolves in 3+1 dimensions accord- 
ing to the conservation laws for energy and momentum, 
i.e. relativistic hydrodynamics. However, the "fluid" of 
quarks interacts locally with the chiral fields, that is, 
they can exchange energy and momentum. In turn, the 
(long wavelength modes of the) fields obey the classi- 
cal equations of motion which follow from the underly- 
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ing Lagrangian in the presence of the quarks and anti- 
quarks. Similar models for the dynamics of quarks cou- 
pled to chiral fields were considered in the past. In [17], 
a background of freely streaming quarks was assumed, 
and the classical evolution of the chiral fields was dis- 
cussed. More realistic dynamical descriptions for the 
quark medium followed shortly, treating them as either 
a relativistic fluid [18], as also envisaged here, or within 
classical Vlasov transport theory [19,20]. Those studies 
focused on a second-order chiral phase transition, or, in 
the presence of explicit symmetry breaking, on a smooth 
crossover. However, it turns out that one can also ad- 
dress first-order chiral phase transitions within the very 
same model, at least in a phenomenological fashion, by 
chosing larger quark- field coupling g [21,22] (see below). 
Integrating out the quarks then leads to an effective po- 
tential exhibiting two degenerate states around T c . 

Here, wc extend the previous work mentioned above, 
and at the same time shift our focus somewhat. Namely, 
the early studies were mainly concerned with the dy- 
namical evolution of the long-wavelength chiral fields, 
and of classical pion production; that is, they mainly 
adressed issues related to the possiblity of forming "Do- 
mains of Disoriented Chiral Condensates" (DCC), as sug- 
gested by Rajagopal and Wilczek, and others [23], see 
also [13,18,19,21]. Our present work puts more empha- 
sis on the dynamics of the heat-bath of quarks, rather 
than on that of the soft modes of the chiral field. Wc 
shall point out qualitative changes in the classical energy- 
momentum flow of the "fluid" of quarks in the proximity 
of a chiral critical point, rather than look for "rare phe- 
nomena" like DCC formation. 

Moreover, refs. [17-21] all employed the mean field ap- 
proximation for the chiral fields. Field fluctuations at the 
phase transition were not considered. As an example, for 
the first-order phase transition discussed in [21] dynami- 
cal bubble nucleation ( "boiling" ) could not be described, 
as it requires large coherent thermal field fluctuations 
from the symmetry restored phase, over the barrier and 
into the symmetry broken phase (see e.g. [24] for results 
of such dynamical simulations, and [22] for a computation 
of bubble nucleation rates from the linear sigma model). 
Thus, the main improvement here is that we do include a 
dynamical treatment of field fluctuations in the vicinity 
of a critical point, and their influence on the dynamical 
evolution of the quark fluid. 

On the technical side, going beyond the mean field 
approximation requires us to introduce appropriate sub- 
tractions in all thermodynamical functions, as explained 
in appendix A. Moreover, the coupled system of non- 
linear partial differential equations has to be solved 
numerically in 3+1 dimensions, without imposing any 
space-time symmetry assumptions (while [17-19,21] all 
simplified the solution greatly by assuming special sym- 
metries which essentially reduced the problem to 0+1 
or 1+1 space-time dimensions). That is because fluctua- 
tions break any space-time symmetry that may be obeyed 
by the mean field, as for example spherical symmetry or 



symmetry under Lorentz boosts in a particular direction. 

As mentioned in the introduction, physically the chi- 
ral critical point is expected to occur for some specific 
values of temperature T and baryon-chemical potential 
fiB- To simplify the problem and its numerical solution, 
however, here we rather choose to consider only locally 
baryon symmetric matter, i.e. equal numbers of quarks 
and anti-quarks. Instead, we can "shift" the critical point 
by varying the quark-field coupling constant g, that is, 
by increasing or decreasing the vacuum mass of the con- 
stituent quarks. As explained in more detail below, large 
values for g result in a first-order phase transition, while 
small g leads to a crossover. For the observables dis- 
cussed in section III, the qualitative difference between 
the two realizations of a chiral critical point mentioned 
above should not matter much 1 . On the other hand, our 
simplified treatment disables us from studying fluctua- 
tions of net baryon charge [26]. 

In section II A we discuss the effective potential "seen" 
by the long wavelength modes of the chiral fields in the 
presence of a heat bath of quarks and anti-quarks. Then, 
in II B we present our model for the non-equilibrium dy- 
namical treatment of field and fluid evolutions. Some 
numerical algorithms and details are mentioned briefly 
in section II C. In section III we present our results and 
end with a summary and an outlook in section IV. 

A. Effective Potential 

As an effective theory of the chiral symmetry breaking 
dynamics, we assume the linear cr-model coupled to two 
flavors of quarks [9] : 

C = q [ij^dp - g{a + 75T • 7?)] q 

+ ^(d tl ad^a + d tl 7fd^7f)-U(a,7f) . (2.1) 

The potential, which exhibits both spontaneously and 
explicitly broken chiral symmetry, is 

U(a,Tr) = ^(a 2 +n 2 -v 2 ) 2 -h q a-U . (2.2) 

Here q is the constituent-quark field q = (u,d). The 
scalar field a and the pseudoscalar field if — (771,772,773) 
together form a chiral field (f> a — (<j,tt). The pa- 
rameters of the Lagrangian arc chosen such that chiral 
SUl(2) ® SUr(2) symmetry is spontaneously broken in 
the vacuum. The vacuum expectation values of the con- 
densates are (a) = f n and (77) = 0, where / w = 93 MeV is 



1 However, we note that there are indeed differences on the 
quantitative level. For example, for a first-order phase tran- 
sition in baryon dense matter the isentropic speed of sound 
does not vanish in general at T c [16,25], as it does for zero net 
baryon charge, pis = 0. 
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the pion decay constant. The explicit symmetry breaking 
term is due to the non-zero pion mass and is determined 
by the PCAC relation, which gives h q — f n rn^, where 
= 138 MeV. This leads to v 2 = f 2 - ml/X 2 . The 
value of A 2 = 20 leads to a cr-mass, m 2 = 2A 2 f 2 +m 2 , ap- 
proximately equal to 600 MeV. In mean field theory, the 
purely bosonic part of this Lagrangian exhibits a second- 
order phase transition [8] if the explicit symmetry break- 
ing term, h q , is dropped. For h q ^ 0, the transition be- 
comes a smooth crossover from the phase with restored 
symmetry to that of broken symmetry. The normaliza- 
tion constant Uq is chosen such that the potential energy 
vanishes in the ground state, that is 



Here, d„ 



rr - m * 



r2 1 



(2.3) 



For g > 0, the finite-temperature one-loop effective 
potential also includes a contribution from the quarks. 
In our phenomenological approach, we shall consider the 
quarks as a heat bath in (local) thermal equilibrium. 
Thus, it is possible to integrate them out to obtain the 
effective potential for the chiral fields in the presence of 
that bath of quarks. Consider a system of quarks and 
antiquarks in thermodynamical equilibrium at tempera- 
ture T and in a volume V. The grand canonical partition 
function reads 



Z 



VqVqVaV'K exp 



l/T 



d(it) 



d 3 xC 



• (2-4) 



(Anti-) Periodic boundary conditions for the (fermion) 
boson fields are implied. In mean-field approximation 
the chiral fields in the Lagrangian are replaced by their 
expectation values, which we denote by a and ir. Then, 
up to an overall normalization factor, 



Z = Mi 



u 



jvqVq exp 



d 3 x 



Q [il^dp ~ giv + il5T ■ 7f)] q , 



= Mu det p { [prf - g{a + i lb f • if)] /t} , (2.5) 



where 



Mu = exp 



VU{a, t?) 



(2.6) 



Taking the logarithm of Z, the determinant of the Dirac 
operator can be evaluated in the standard fashion [27], 
and we finally obtain the grand canonical potential 



--\ogZ = U + V q (T) , 



(2.7) 



where 



V q {T) = -d q J j^{E + T log 



1 + e 



-E/T 



]} ■ (2-8) 



24 denotes the color-spin-isospin-baryon 



charge degeneracy of the quarks. For our purposes, the 
zero-temperature contribution to V q , i.e. the first term in 
the integral in (2.8), can be absorbed into U via a stan- 
dard renormalization of the bare parameters A 2 and v 2 . 
The logarithmic dependence on the renormalization scale 
is ignored in the following. 

Adding the T — and the finite-T contributions de- 
fines our effective potential V e g: 



K ff (0 a ,T) = C/(0 a ) 



d q T 



d 3 p 
(2tt)3 



log (l + e- E ' T ) 



(2.9) 



V e B depends on the order parameter field through the ef- 



fective mass of the quarks, m q = g\ 
which enters into the expression for the energy, E 

\fp 2 + g 2 <t> 2 - 

In principle, one should also integrate out short wave- 
length fluctuations of <f), which would lead to an addi- 
tional contribution to the effective potential (2.9), see for 
example [13,28]. For the present phenomenological anal- 
ysis, however, the expression (2.9) is already sufficient 
in that it exhibits a critical endpoint for some particular 
value of the coupling constant g. Since we do not ex- 
pect the simple model (2.1) to be quantitatively reliable 
anyway, it is not unreasonable to employ the effective po- 
tential (2.9) for a study of qualitative effects near a chiral 
critical point. 




0.4 0.5 

FIG. 1. The one-loop finite-T effective potential as a func- 
tion of the scalar field a at if = 0. The quark-field coupling 
constant g is being varied. The curves are labeled by the 
value for g and by the temperature in MeV in parantheses. 



The field self-coupling is chosen to be A 



20. 



We now turn to a discussion of the shape of the effec- 
tive potential, cf. Fig. 1. For sufficiently small g one still 
finds the above-mentioned smooth transition between the 
two phases. At larger coupling to the quarks, however, 
the effective potential exhibits a first-order phase tran- 
sition [21]. Along the line of first-order transitions, for 
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temperatures near the critical temperature, V e g displays 
a local minimum a = a\ (T) ~ which is separated by a 
barrier from another local minimum at a — 02 (T) > 0. 
(There is another local minimum for negative a which is 
of higher energy and does not concern us.) These two 
minima are degenerate at T = T c . For example, g = 5.5 
leads to a critical temperature of T c ~ 123.3 MeV. Low- 
ering the value of g leads to a smaller barrier between 
the two degenerate states. Also, o\ approaches 02, i.e., 
the phase transition weakens, and moreover the spinodal 
temperature approaches T c [22]. At g c ~ 3.7, finally, 
the barrier disappears, and so the latent heat vanishes. 
This is the second-order critical point, where the poten- 
tial about the minimum is flat. 

A different possibility of making the cr-meson much 
lighter at T c than at T = is to reduce the self-coupling 
of the chiral fields [29], A 2 , rather than that to the quarks. 
For example, one may choose A 2 ~ 2.2 with the pion de- 
cay constant f n and vacuum mass fixed, such that 
v 2 = 0. However, within the present model this also re- 
duces T c significantly, to less than 100 MeV [22]. Such 
low T c appear to be excluded by present lattice QCD 
results [30] , and moreover would generate too small ther- 
mal fluctuations in the heat bath. Therefore, we keep 
A 2 = 20 fixed throughout the manuscript. 



B. Coupled Dynamics of Fields and Fluid 



g(qq) = 
g{qi5Tq) = 



5 (n q ) _ 5 (v cS - u) 

5a 5a 
5 (H q ) _ 5 (V eS - U) 



Sir 



5tt 



(2.13) 



Applying this to the right hand side of equation (2.9) 
yields the expressions for the scalar density and the pseu- 
doscalar density as given in (2.11). 

As already mentioned above, we shall assume that the 
quarks constitute a heat bath in local thermal equilib- 
rium. Thus, their dynamical evolution is determined by 
the local conservation laws for energy and momentum in 
relativistic hydrodynamics. For simplicity, we shall fur- 
ther assume that the stress-energy tensor of the quark 
fluid is of the "perfect fluid" form (corrections could in 
principle be taken into account in the future along the 
lines discussed in [31]), 



Here, 



T"" = {e + p)uV - pg^ 



T» v u v 



^Ju a T°PT pa u 



(2.14) 



(2.15) 



is the local four-velocity of the fluid. g^ v = 
diag(l, — 1, — 1, — 1) is our metric tensor, and so the line 
element is ds 2 = dt 2 — dx 2 . The time t is measured in 
the global restframe. Furthermore, 



The classical equations of motion for the chiral fields 
are 



d^a +^- = -g (qq) = -gp s 
da 



5U 

d ^ +^ = -9 i^rq) = -gp P s, (2.10) 



where 



Ps = (qq) =9<7d q J ?J^3^/(P) » 



(2tt) 3 £" 
f d 3 p 1 

P P s = (qibfq) =gnd q J J^TagfiP) ( 2 - n ) 

are the scalar and pseudoscalar densities generated by the 
heat bath of quarks and anti-quarks, respectively. The 
distribution of the quarks and anti-quarks in momentum 
space is given by the Fermi-Dirac distribution. 

The form (2.10) for the coupling to the heat bath of 
quarks can be derived from the Lagrange density (2.1) 
in mean-field approximation. The energy density of the 
quarks is given by 



e(4>,T) = (H q ) , 

p(cb,T) = -V cB (cb,T) + U(cf, a ) 



(2.16) 



denote the energy density and the pressure of the quarks 
at temperature T. Note that we do not assume that 
the chiral fields are in equilibrium with the heat bath 
of quarks. Hence, both e and p depend explicitly on 
(f) = v 'a 2 + 7? 2 . Given an initial condition on some space- 
like surface, at any other causally connected space- 
time point is determined by 



S v 



(2.17) 



In the absence of interactions between the chiral fields 
and the quarks, the source term 5"^ vanishes, and the 
energy and the momentum of the quarks, 



(E,P)= / da^ 



are conserved: 



u^E = u^P 1 = 



(2.18) 



(2.19) 



(Hq) = (qi-h) + 9 (qoq) + g (qi^f ■ nq) 



(2.12) 



The finite-T contribution to the equation of motion for 
the a, 7? fields is obtained from the variation of the effec- 
tive potential with respect to a or 7? , respectively: 



With interactions turned on, this is obviously not true 
any longer. Rather, the total energy and momentum of 
fluid plus fields is the conserved quantity: 



(2.20) 
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The stress-energy tensor for the fields can be computed 
from the Lagrange density in the standard fashion. The 
effective mass of the quarks, m q = g\/^2 a is al- 
ready accounted for in the EoS for the quark fluid, see 
eqs. (2.16). Thus, T£" is the stress-energy tensor of the 
chiral fields alone: 

£* = £|(8 1 A)(3"*.)-tf ■ < 2 ' 21 ) 

a 

Its divergence is given by 

= g Ps d v a + gp ps ■ d"n . (2.22) 

In the second step we made use of the equations of mo- 
tion (2.10). This is the source term in the continu- 
ity equation (2.17) for the stress energy tensor of the 
quark fluid. A different derivation based on moments of 
the classical Vlasov equation for the quarks was given 
in [17,18], see also [19]. 

We emphasize again that we employ eqs. (2.10) not 
only to propagate the mean field through the transition 
but fluctuations as well. The initial condition includes 
some generic "primordial" spectrum of fluctuations, see 
section III A, which then evolve in the effective potential 
generated by the matter fields, i.e. the quarks. Near the 
critical point, those fluctuations have small effective mass 
and "spread out" to probe the flat effective potential. 
Since all field modes are coupled, effectively the fluctu- 
ations act as a noise term in the equation of motion for 
the mean field, similar to the familiar Langevin dynam- 
ics [32]; near the phase transition, however, the noise is 
neither Gaussian (the effective potential is not parabolic) 
nor Markovian (zero correlation length in time) nor white 
(zero correlation length in space). Rather, correlation 
lengths and n-point functions of the "noise" are governed 
by the dynamics of the fluctuations in the effective po- 
tential generated by the fluid of quarks, including also 
effects from the finite size and the relativistic expansion 
of the system. 

C. Numerics and Technical Details 

We briefly describe how we solved the coupled sys- 
tem (2.10,2.17) of partial differential equations in 3+1 
space-time dimensions. 

We follow the evolution on t = const hypersurfaces, 
and in a fixed spatial cube of volume L 3 . We discretize 
three-space in that cube by introducing a 160 3 grid with a 
spacing of Ax = 0.2 fm (thus, L = 160 x 0.2 fm=32 fm). 
On that grid, we solve the hyperbolic continuity equa- 
tions of fluid dynamics (2.17) using the so-called phoeni- 
cal SHASTA flux-corrected transport algorithm with 



simplified source treatment. It is described and tested 
in detail in Ref. [33], and we refrain from a discussion 
here. The time step was chosen as At = 0.4Ax, as ap- 
propriate for the SHASTA [33,34]. We performed each 
time step in the standard fashion with S u = 0, then sub- 
tracted the sources S v from the energy and momentum 
density in the calculational frame (i.e. the global rest 
frame of the fluid). Finally, from T 00 , T oi and the EoS 
p = p(e, c/) a ) we solved for the velocity of the local rest- 
frame (LRF) of each cell, and for the energy density e of 
the fluid in the LRF. Such a treatment of sources in the 
continuity equations for energy and momentum proved 
to work well in relativistic multi-fluid dynamics [16,34], 
where one encounters a similar equation due to interac- 
tions between various fluids. The boundary conditions 
at the edge of the computational grid are such that the 
fluid simply streams out when reaching the boundary. 
This can be monitored by checking the conservation of 
the total energy as a function of time. 

Regarding the fields, we solve the classical equations 
of motion using a staggered leap-frog algorithm with 
second-order accuracy in time, see for example Ref. [35]. 
The non- linear wave-equations (2.10) are split into two 
coupled first-order equations (in time) by considering 
separately the field <p a and its canonically conjugate field 
d(j)a/dt. For this algorithm, the time step for propa- 
gation of the fields has to be chosen smaller than that 
for propagation of the fluid. In practice, we found that 
At = 0.1 Ax gave reasonable accuracy. As the time 
step for propagation of the fluid was chosen to be four 
times larger, the fields where propagated for four con- 
secutive steps, with the fluid kept "frozen", followed by 
one single step for the fluid. Spatial gradients of the 
fields where set to zero at the edges of the computa- 
tional grid, rather than employing periodic boundary 
conditions. The spatial derivatives in the equations of 
motion for the fields (2.10) where discretized to second 
order accuracy in the grid spacing to match the second 
order accuracy in time. We employed the same small grid 
spacing Ax = 0.3 fm as for the fluid. Although the mean 
fields vary on a larger scale, it is important to allow for 
non-linear amplification of harder fluctuations which can 
couple to the soft modes and affect the dynamical relax- 
ation to the vacuum state with broken symmetry, see for 
example [36] and [37]. 

For the fluid, one needs to introduce nine 3d spatial 
grids 2 , for e, p, T° M , and u. In addition, one needs two 
more 3d grids for each field component, namely for </> a 
and d(j)a/dt. To save computational resources, we have 



2 At finite baryon density, one would need another grid for 
Pb- Also, inverting the Lorentz transformation to obtain the 
LRF densities e and ps from T 00 , T (H and the net baryon 
current Jg = pbu 11 in each time step would require a three- 
dimensional root search for the function p = p(e, Pb,4>)- 
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therefore decided to neglect the pion field altogether, i.e., 
to set it to zero everywhere in the forward light cone. 
Thus, only the cr-field is considered in the actual com- 
putations described below. However, this represents a 
minor restriction only, since we focus here on the bulk 
evolution of the quark fluid rather than on fluctuation 
observables or coherent pion production from the decay 
of a classical pion field. (Coherent pion production is 
known to contribute a small fraction of the total pion 
yield only [18,19,21,23,38]. Incoherent particle produc- 
tion from the decay of the fluid by far dominates, if no 
kinematic cuts are applied.) 



III. RESULTS 

Having formulated our model, we now proceed to show 
some specific examples of numerical solutions. In partic- 
ular, we would like to examine whether the dynamical 
evolution changes as one crosses the chiral critical point. 



A. Initial Condition 

As an example, we employ the following initial con- 
ditions for the results described below. Of course, one 
could employ more refined initial conditions and "tune" 
them such as to reproduce various aspects of the final 
state, which in principle could be compared to experi- 
mental data. At present, our more modest goal is to 
illustrate qualitative effects originating from the phase 
transition. The initial conditions are meant to provide a 
semi- realistic parameterization of the hot fireball created 
in a high-energy heavy-ion collision. 

At time t — 0, the distribution of energy density for 
the quarks is taken to be uniform in the z-direction (with 
length 21 z = 12 fm) and ellipsoidal in the x — y-plane: 



e(t = 0,x) 



s oq : x 2 b 2 + y 2 a 2 < (ab) 2 A z < l z 
: x 2 b 2 + y 2 a 2 > (ab) 2 V z > l z 



(3-1) 



where e oq denotes the equilibrium value of the energy 
density taken at a temperature of Tj « 160 McV. In our 

calculation we choose a = ta — b/2 and b = \J r 2 A — 6 2 /4, 

where b denotes the impact parameter of two nuclei with 
radius ta (below, we choose ta = 6.5 fm and b = 6 fm). 
Thus, the ellipsoidal shape resembles the almond shaped 
overlap of two colliding nuclei. The 9-function distribu- 
tion of the initial energy density is evidently somewhat 
unrealistic; a smoother distribution with non-zero surface 
thickness would be more realistic and perhaps affect the 
results somewhat. However, as already mentioned above, 
here we do not aim at quantitative fits to experimental 
data but at illustrating qualitative effects related to the 
shape of the effective potential in the transition region. 



The collective longitudinal velocity of the fluid of 
quarks is assumed to rise linearly with z: v z (t = 0,x) oc 
z/l z ■ v maxi where v max = 0.2 . The transverse compo- 
nents of v are set to zero at t = 0. 

Our initial conditions for the chiral fields are 



a(t = 0,x) = Sa(x) + U + 

( — fir + <7eq) ' 



• 1 + exp 

7?(t = 0, x) = Sn(x) , 



1 + exp 



R 



where f = \J x 2 + y 2 , 



abf 



R = 



\Jb 2 x 2 + a 2 y 2 



r^0 
f = 



(3.2) 



(3.3) 



and a = 0.3 fm is the surface thickness of this Woods- 



Saxon like distribution. Here a P 



is the value of the 



a field corresponding to e oq . Thus, the chiral condensate 
nearly vanishes at the center, where the energy density 
of the quarks is large, and then quickly interpolates to 
f v where the matter density is low. 

5a(x) and 5tt(x) represent the initial random fluctu- 
ations of the fields. Our focus at this stage is on how 
those "primordial" fluctuations evolve through the phase 
transition (or crossover) and how they affect the hydro- 
dynamic expansion of the thermalized matter fields (the 
fluid). Thus, for a first qualitative analysis we do not rely 
on additional physics input 3 for the primordial fluctua- 
tions but rather chose a generic Gaussian distribution, 



P[Stj> a ] ex exp (-501/2 (5<g)) 



(3.4) 



with the variance (Scf) 2 ) left as a free parameter. (Here, 
no summation over the index o (internal space) is car- 
ried out.) We performed simulations with yj (<50g) = 
\J (<5c 2 ) = v/3. Such fluctuations are large enough to 
probe the barrier of the effective potential in case of a 
first order transition, or the flat region close to the criti- 
cal endpoint (Fig. 1). On the other hand, they are suffi- 
ciently small to allow for a one-loop subtraction of their 
contribution to the effective potential, as described in 
appendix A. 



3 For example, one might assume thermalized primordial 
fluctuations, in which case their distribution depends on the 
effective potential at the temperature T;. We have checked, 
however, that at high temperature V e g (<j>, T{) looks rather sim- 
ilar for all values of g considered here, regardless of whether 
later on the evolution proceeds through the crossover, the 
critical endpoint or the first-order phase transition regimes. 
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The time derivatives of the fields where set to zero 
at t = 0, with no fluctuations. As already mentioned in 
section II C, the actual numerical computations described 
below were performed with 5w = 0, i.e. (84>\ 2 3) = 0. 

One must further take into account that the field 
fluctuations are correlated over some space-like distance 
£ ~ 1 fm. This is a physical scale which is present in the 
initial conditions; if one simply picks random fluctuations 
at each point of the grid, the correlation length will in- 
stead be given by the artificial numerical discretization 
scale Ax. 

In practice, a useful and simple procedure to imple- 
ment physical correlations in the initial condition is as 
follows. First, at each point of the grid one samples the 
distribution (3.4) at random. Then, one smoothly sweeps 
a " sliding window" of linear size £ = nAx over the grid 
and averages the fields: 




x 4> a (x + iAxei + j Ax<32 + kAxe^) . (3.5) 

i : j,k— 0. .71— 1 

Here, e\, &2, S3 define a global cartesian orthonormal ba- 
sis, and x = le*i + IS2 + le*3, lei + le*2 + 2e3, • • •, is the 
sequence of grid points. Clearly, the above averaging pro- 
cedure "cools" the fluctuations, in that (Stfi'^) ^ {$4>V)- -"- n 
particular, in the continuum limit Ax — ► with £ = nAx 
held constant, one of course finds that (5<f>'a) ~ > for the 
distribution (3.4). Therefore, in a final sweep over the 
grid one has to rescale the fields at each grid point by 




This procedure leads to an initial field configuration that 
exhibits both the proper physical correlation length and 
the desired fluctuations. Moreover, it prevents the initial 
energy density from field gradients to grow like 1 /Ax 2 . 

B. Numerical Solution 

In the following we consider both a first order phase 
transition corresponding to g — 5.5, as well as the critical 
point at g = 3.7. 

Figs. 2 and 3 depict the time evolution of the a field 
along the x-axis and the y-axis, respectively. At t = 
the field within the hot region has small amplitude, 
corresponding to the chiral symmetry restored phase. 
That region is surrounded by the physical vacuum with 
(a) = f n = 93 MeV. 



9 = 3-7 I 




x [fm] 



FIG. 2. Space-time evolution of the chiral field along the 
rs-axis at y — z = 0. Dark regions correspond to small field 
amplitudes (symmetric phase), light regions to large ampli- 
tudes (broken phase). The scale on the right specifies the 
field amplitude in MeV. 

For the first order phase transition (g = 5.5) a bar- 
rier separates the two degenerate minima of the effective 
potential (Fig. 1) at T c . Figs. 2 and 3 show that this 
barrier leads to a rather well-defined surface in coordi- 
nate space, separating the vacuum from the symmetric 
phase. For the above-mentioned initial conditions, the 
fluctuations are not strong enough for the field to eas- 
ily overcome the barrier. Nevertheless, one can observe 
dynamical fluctuations into the broken symmetry state, 
e.g. at x ~ 1 fm and t ~ 2 — 4 fm/c, which however col- 
lapse again. Our dynamical results agree with previous 
arguments that nucleation is a slow process on the time 
scale of heavy-ion collisions, and so the Gibbs phase equi- 
librium is not established dynamically [21,22,24,39]. At 
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time t ~ 9 fm the phase transition occurs spontaneously 
"in an instant" , that is, on a space-like surface which 
can clearly be seen in Figs. 2 and 3. A "bubble" of the 
symmetric phase survives at x ~ 3 fm for a long time. 




shows a histogram of the field distribution at the center 
(\/ x 2 + y 2 + z 2 < 2 fm). One observes that the distri- 
bution broadens from t = 4 fm to t = 10 fm, and then 
narrows again at later times after the transition to the 
broken phase occured. Also, notice that the distributions 
at t = 4 fm and t = 14 fm are well described by Gaus- 
sians, i.e. the effective potential is essentially parabolic, 
while for t = 10 fm there are visible deviations from a 
simple Gaussian. 
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FIG. 4. Histogram of the field distribution at the center 
for three different times. Lines represent Gaussian fits with 
standard deviations 25 MeV (t=4 fm), 53 MeV (t=10 fm), 
and 43 MeV (t=14 fm). 



3 4 
y[fm] 

FIG. 3. Space-time evolution of the chiral field along the 
y-axis at x = z — 0. 



The picture is rather different for the transition at the 
critical point, i.e. g = 3.7. Here, the barrier between the 
degenerate minima vanishes and the potential is flat. As 
is evident, there are no clear surfaces separating either 
the vacuum from the center or high-density bubbles (or 
"droplets" ) from their surrounding. Due to the flatness 
of the potential, near the center the field performs large- 
amplitude oscillations (from a ~ to a > / w ) for a long 
time; they extend in space over distances k, 1 — 3 fm 
(e.g. at t w 9, 12, and 14 fm/c in Fig. 2), which is not 
much less than the initial size of the hot region. Figure 4 
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g = 3.7 | 
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FIG. 5. Space-time evolution of the fluid energy density 
along the x-axis at y = z = 0. The scale on the right spec- 
ifies the energy density in units of nuclear matter density 
e « 150 MeV/fm 3 . 

The time evolution of the local rest-frame energy den- 
sity e of the quarks is shown in Figs. 5, 6. Again we see 
large-scale structures for the first order phase transition 
(g = 5.5), while the energy density is rather homogeneous 
on large time and distance scales if the expansion trajec- 
tory goes through the critical endpoint. For the first- 
order transition, quarks can be "trapped" in droplets 
with (7^0 (the minimum of the effective potential where 
the symmetry is restored) because the mass barrier can 
keep them from escaping. The droplet of high-density 
matter at x ~ 3 fm and i ~ 12 — 16 fm/c can easily be 
associated to the region of nearly vanishing chiral scalar 
field from Fig. 2. Eventually, that region must perform 
the transition to the symmetry broken state, either by a 
strong thermal fluctuation or when reaching the spinodal 



point. At the spinodal, the system is as far from local 
thermal equilibrium as it can get, and the "roll-down" of 
the order parameter field to the global minimum of the 
potential can influence the collective expansion of the 
quark fluid. 

9 = 3-7 | 




3 4 
y[fm] 

FIG. 6. Space-time evolution of the energy density along 
the y-axis at x = z = 0. 

Note that the energy density at the center drops more 
rapidly for the first-order transition than near the criti- 
cal endpoint. This has consequences for the build-up of 
azimuthally asymmetric flow, as we shall discuss below. 

Fig. 7 depicts the time evolution of the azimuthal mo- 
mentum anisotropy [40] 



(T X X Tyy) 

(T x: 



T, 



(3.7) 



Hill 



where the averages of the stress-energy tensor of the fluid 







are taken at fixed time: 



0.25 



(Tij){t) = / d 3 xT i:j {t,x) 



(3.8) 



Also, we average over a few initial field configurations, 
which gave similar results for e p (t), though. 

For the above-mentioned initial conditions the energy- 
momentum tensor is symmetric, and so e p = at t = 
(this might be different in more realistic treatments [41]). 
Pressure gradients in x— and y— directions are differ- 
ent, though. Therefore, the acceleration of the fluid is 
stronger in the reaction (x — z) plane than out of plane, 
leading to a nonzero azimuthal asymmetry e p > at 
times t > 0. The asymmetry first grows nearly linearly 
with time but saturates when the asymmetry of the en- 
ergy density and of the pressure gradients becomes small. 
As explained above (Figs. 5, 6), this happens earlier for a 
first-order transition than for trajectories near the chiral 
critical endpoint. This is then reflected in the final value 
of e p . We stress that the more rapid saturation of the 
azimuthal asymmetry in case of a first-order transition 
is not in contradiction to the fact that hot (high-energy 
density) "droplets" survive for rather long times, as seen 
in the figures. Rather, such "droplets" typically turn out 
to be more or less rotationally symmetric, or at least ex- 
hibit deformations which are uncorrelated to the reaction 
plane (the x — z plane in our case). Thus, they tend to 
reduce the average azimuthal asymmetry of the energy- 
momentum tensor. 

For comparison, in Fig. 7 we also show the result for an 
equilibrium first-order phase transition. Here, the equa- 
tions of motion for the chiral fields, eqs. (2.10) are not 
solved but rather the tr-field is required to populate the 
(global) minimum of the effective potential, 



SVcB 

5a 
S 2 V e s 
Sa 2 



= 



> 



(3.9) 
(3.10) 



That is, the chiral field is in equilibrium with the quark- 
antiquark fluid and does not exhibit any explicit space- 
time dependence. At T c , where two degenerate minima 
exist, one performs the usual Maxwell- Gibbs construc- 
tion to determine the fractions of the total volume occu- 
pied by matter in the symmetric and the broken sym- 
metry phases, respectively. Evidently, for the above- 
mentioned initial conditions the equilibrium phase tran- 
sition leads to nearly the same azimuthal asymmetry of 
T IJ as for the cross over. Therefore, it is indeed the non- 
equilibrium real-time dynamics (field fluctuations over 
the free-energy barrier) that is responsible for the ob- 
served reduction of e p in the regime of first-order chiral 
phase transitions. 
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— g = 3.7 

— g = 5.5(eq.) 

— - g = 5.5 
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FIG. 7. Time evolution of the momentum anisotropy for a 
cross over (g — 3.3), a second-order phase transition at the 
critical point (g = 3.7), and a first order phase transition 
(g — 5.5) in or out of equilibrium. 



IV. SUMMARY AND OUTLOOK 

In summary, we have introduced a simple phenomeno- 
logical description of the non-equilibrium real-time dy- 
namics of the chiral phase transition in an expanding 
(relativistic) fluid of quarks. More precisely, we coupled 
the linear sigma model, which describes the dynamics of 
the long-wavelength modes of the chiral order parame- 
ter field, to the hydrodynamical evolution of a system of 
quarks. The chiral field(s) evolve according to the finite- 
temperature effective potential that is generated by in- 
tegrating out the quarks from the Lagrangian; in turn, 
the field(s) determine the effective quark mass (i.e. the 
equation of state of the quark fluid) dynamically. 

The above model exhibits a first order phase transi- 
tion for large g, which is the quark- field coupling con- 
stant. The line of first order transitions ends in a critical 
point when g is lowered, i.e. the transition turns into a 
cross over for smaller couplings. Thus, by varying g one 
can qualitatively compare the hydrodynamic expansion 
pattern of the quark fluid for dynamical trajectories that 
cross the line of first order transitions to that obtained 
in the cross over regime. 

We have obtained numerical solutions in 3+1 space- 
time dimensions, using simple initial conditions that 
might be appropriate for relativistic heavy-ion collisions. 
The hydrodynamical expansion pattern clearly depends 
on the structure of the effective potential. For trajecto- 
ries in the cross over regime or near the critical endpoint 
the overall bulk dynamics is found to be rather "smooth" , 
in that the space-time distribution of the energy density 
of the fluid is not affected very much by the fluctuations 
of the order parameter field. In the absence of a latent 
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heat, the energy density can not jump much between re- 
gions where the field amplitude is different. 

In contrast, if the effective potential exhibits a bar- 
rier between the symmetry restored and broken phases, 
respectively, we do see that large-scale structures are 
formed dynamically, e.g. "droplets" of the symmetric 
phase may survive for rather long times before becoming 
mechanically unstable (at the spinodal). In that sense, 
the overall time scale is longer for trajectories crossing the 
line of first order transitions. Nevertheless, typically such 
structures are not correlated to the reaction plane; thus, 
the direct correspondence of spatial anisotropies in the 
initial condition to momentum-space anisotropies in the 
final state predicted by equilibrium hydrodynamics (that 
is, when the phase transition is not treated dynamically 
but modelled by a Maxwell- Gibbs construction) is weak- 
ened. For example, we find much smaller momentum- 
space anisotropy for a dynamical first-order transition 
than for a trajectory through the chiral critical endpoint 
(for the same initial condition). This could be a very 
useful prediction with regard to the experimental search 
for the chiral critical endpoint of QCD in heavy-ion colli- 
sions at the BNL-AGS, the CERN-SPS and the envisaged 
new GSI heavy-ion accelerator. Until now, experiments 
focused on fluctuation observables, but inclusive observ- 
ables usually are much easier to analyze accurately. 

In the future, we intend to scrutinize other inclusive ob- 
servables as to their sensitivity to non-equilibrium effects 
from phase transitions. Of course, there is also plenty of 
room to improve on the model in order to obtain more 
quantitative predictions. The present paper represents a 
first step towards an actual real-time description of the 
chiral phase transition on either side of the critical end- 
point in expanding relativistic fluids with realistic 3+ Id 
geometries. 



Note added: After this manuscript was submitted for 
publication the NA49 collaboration published the elliptic 
flow at £i ab = 40,4 GeV [42]. From Fig. 24 of that publi- 
cation, the dependence of v 2 on log ^/s is approximately 
linear. However, the "natural" scale for v 2 is set by (p t ), 
not log-^s, as pointed out by Snellings [43]. Indeed, at 
high energies the differential v 2 (p t ) of charged hadrons 
is approximately proportional to p t , such that the aver- 
aged v 2 oc (pt)- In fact, for mid-central collisions v 2 in- 
creases from w 3% at top SPS energy (y/s — 18A GeV) to 
w 4.5% at RHIC energy (^/i = 130,4 GeV). When scaled 
by the average transverse momentum, though, the ellip- 
tic flow in that energy regime is nearly constant [43] . To 
scrutinze deviations from the "natural" scaling v 2 oc (pt), 
we plot the excitation function of V2/(pt) in Fig. 8. 
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FIG. 8. Excitation function of V2 / (pt) of negatively charged 
particles in mid-central collisions from top AGS to RHIC en- 
ergy. Data for V2 are taken from Fig. 24 in [42] and (p t ) 
from [44-47]. 

One observes that, as already mentioned above, the 
data is compatible with no energy (or (p t )) dependence 
above top SPS energy. Clearly, there is a systematic drop 
of v 2 relative to (p t ) towards lower energies. For instance, 
at £ lab = 40A GeV corresponding to (p t ) ss 350 MeV, 
v 2 /(pt) is lower by about two standard deviations than 
at higher energies. Qualitatively, this interesting behav- 
ior is similar to the reduction of the azimuthal momen- 
tum asymmetry, predicted above, caused by crossing the 
second order critical point into the regime of first order 
phase transitions. Additional studies of "conventional" 
non-equilibrium effects unrelated to a phase transition 
are certainly required, however, before firm conclusions 
can be drawn. 
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APPENDIX A: SUBTRACTING INITIAL 
FLUCTUATIONS 

In section II A we discussed our phenomcnological 
ansatz for the effective potential for the long wavelength 
modes of the chiral fields, as generated by the heat bath 
of quarks. Formally, it is obtained from the Lagrangian 
by integrating out the quarks to one loop. 

Our main objective here is to study dynamically fluctu- 
ations of the chiral order parameter (or effects generated 
by those fluctuations) in the vicinity of the chiral critical 
point as the system makes the transition to broken chiral 
symmetry at low temperature. Thus, we have to allow 
for primordial fluctuations of the chiral fields, also. How- 
ever, those fluctuations of the fields at time t — will of 
course also contribute to the effective potential and "dis- 
tort" its shape. In order to restore our original choice for 
V e ff from eq. (2.9), and thus ensure the correct dynam- 
ics for the long-wavelength modes, we have to introduce 
appropriate subtractions. 

The procedure is as follows. The scalar and pseu- 
doscalar densities are given by eqs. (2.11). They depend 
explicitly on the value of the fields (fi a , which we formally 
separate into short and long wavelengths: 



b a (x) = ((pa) + 84> a (x) . 



(Al) 



Here, (•) denotes a spatial average over a volume large 
enough for the fluctuations to average out: 



(Ha) = . 



(A2) 



The linear dimension of that volume will be given by the 
wavelength of the soft modes of interest. 

We now substitute (Al) into eqs. (2.11) and perform 
an expansion up to second order in 5(j) a (x). We then 
perform the averaging over the fluctuations with the dis- 
tribution (3.4) and obtain 



(ps 



Ps((<t>)) 

E 

a 

+ 



I) 



d 3 k 

(2^)3 

9 2 f(k) 
E 2 T 



| + /(fc)e £/T ) 



9 



E 3 T 2 
T 



\f(k) ( 2/ 2 (fc)e 2 ^ - f(k)e E ' T 



+ 3-/(^ + 3^ 



- 2 {5a 2 ) 



2 \ 9 2 f(k) 



E 2 T \E 



+ f(k)e 



E/T 



(A3) 



Here, (S(j) 2 ) = ^2 a (<50 a <50a) is the variance of the fluctua- 
tions in the initial condition, i.e. at the initial time t = 0, 
summed over internal quantum numbers. We made use of 
the fact that the fluctuations (3.4) are diagonal in inter- 
nal space, i.e. (5(pa5<pb) — if a ^ b. The second term is 



the additional contribution seen by the long wavelength 
modes (<fi a ), which is due to the fluctuations. To restore 
the original effective potential, we have to subtract that 
term, i.e. redefine the scalar density as 



Ps(<P,T) = gnd q 



- ^9°d q 



d 3 k 1 

d 3 fc 
(2^p 



f(k) 



g 2 /(fc) 
E 2 T 



f(k)e 



E/T 



+ 



E 3 T 2 



%f(k) (2/W^-/(fc)e^ 



3|/(fc)e £/T 



T 2 \ 
3 2*0 



- 2 (^ 2 )^P(| + /( fc ) e£/T )} (A4) 

This expression has to be substituted for p s on the right- 
hand-side of the equation of motion (2.10). The sub- 
tracted term cancels the "distortion" of the scalar density 
caused by using the local values <r(x), tt(x) for the fields, 
rather than their long- wavelength components (a(x)), 
(n(x)). 

Along the same lines one derives the following expres- 
sions for the fluctuation-subtracted pseudo-scalar den- 
sity, and for the pressure of the quarks: 

f d 3 k 1 
PvM'T) = 9^d q I 77^^;f(k) 



1 



2 
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(2tt) 3 E 
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(2/ 2 (fc) 



2E/T 
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9 2 .f(k) 
E 2 T 
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Izm 1 9 2m 
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Ml 

E 3 T 



f(k)[T 



E 

E-Ef(k)} 



(A5) 
-Sp/Scr and 



One can verify that the identities gp s 
gPps = —Sp/Sir are satisfied, as it should be. 

At fixed values for the fields, the energy density of the 
quarks at a temperature T is given by 



e{<j>,T) 



, dp(cj>,T) 
dT 



■p(4>,T) . 



(A6) 
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Using the expression for the fluctuation-subtracted pres- 
sure given above one obtains 



'g 2 f(k) 



ET 



T - Ef(k)e E ' T 



„4 j.2 



£f(k)[T + E-Ef(k)] 



+ 



E 3 T 



(A7) 



The source term (2.20) changes due to the fluctua- 
tions and one has to use the modified scalar (A4) and 
pseudoscalar (A5) densities, respectively. 

To second order in the fluctuations, the self-interaction 
of the chiral fields is renormalized as 



U(4>a) = y {<? 2 + tt 2 - v 2 ) 2 - h q a - U 

- jE(^) A2 ( 2 ^ + - 2 + 7?2 -- 2 ) 



(A8) 



The above expressions for U((f> a ), p s , P P s, e, and p are 
to be used in the equations of motion for the chiral 
fields (2.10), in the stress-energy tensor of the quark 
fluid (2.14), and in the source term for its divergence 
S v (2.20,2.22). We point out that we subtract those 
quantities for the contribution from initial fluctuations 
of <f> only up to second order in S(f>. We can therefore 
not employ initial conditions with very large local fluc- 
tuations about the mean field. 
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